Dynamical control of the constraints growth in free evolutions of Einstein's equations. 
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I present a new, simple method to dynamically control the growth of the discretized constraints 
during a free evolution of Einstein's equations. During an evolution, any given family of formula- 
tions is adjusted off the constraints surface in a way such that, for any chosen numerical method 
and arbitrary but fixed resolution, the constraints growth can be minimized with respect to the free- 
dom allowed by the formulation. In particular, provided there is enough freedom, the discretized 
constraints can be maintained close to its initial truncation value for all times, or decay from it. 

No a priori knowledge of the solution is needed, and the method can be applied to any formulation 
of Einstein's equations without affecting hyperbolicity. This method is independent of the numerical 
algorithm and accounts for constraint violating modes introduced both by continuum instabilities 
of the formulation and by the numerical method. 



I. MOTIVATION AND OVERVIEW 

The Einstein equations are typically solved in what 
is called free or unconstrained evolutions. The equations 
are split into a set of evolution equations and another one 
of constraints. The standard procedure is to use initial 
data that satisfy the constraints and later solve only for 
the set of evolution equations. The logic is that by virtue 
of Bianchi identities the constraints should be satisfied at 
all times if they do so initially and the evolution equations 
are solved. This is actually true only in the absence of 
boundaries, or in the future domain of dependence of the 
initial data; in the presence of boundaries appropriate 
boundary conditions have to given in such a way that the 
associated initial-boundary value problem is well posed, 
a highly non-trivial issue on its own (see for related 
work) . 

However, problems arise in free evolutions even ne- 
glecting the presence of boundaries. Although the con- 
straints should be exactly satisfied at the continuum 
level, in typical, fixed resolution, numerical simulations 
of strong fields the constraints quickly grow and eventu- 
ally the code crashes. For a numerically stable scheme 
(throughout this paper, the term numerical stability is 
used as equivalent to convergence, in the sense of Lax's 
theorem this growth should go away with resolu- 

tion, but very high resolutions are usually needed if one 
wants to keep these errors under control by just adding 
more points to the simulation. Therefore this procedure 
is not practical and it might not be even feasible. 

Sometimes it happens there are no growing modes 
at the continuum, but the discretization one is using, 
even if numerically stable, can introduce errors that grow 
quickly in time. Examples of this potential source of in- 
stabilities are shown in Ref. ; it is there also explained 
how to rearrange the semidiscrete equations in order to 
prevent this when there is a "conserved" quantity at the 
continuum (for example, in the case of wave propagation 
on a stationary spacetime the physical energy would play 
this role). 

The picture that has emerged in other cases is that 



rapidly growing discrete constraint violating modes are 
numerical excitations of unstable continuum modes. In 
other words, solutions of Einstein's evolution equations 
that are initially slightly off the constraint surface, but 
then deviate from it very quickly. The standard approach 
is then to seek formulations that are as stable as pos- 
sible under continuum constraint violations. A trivial 
observation here is that when rapidly growing discrete 
constraint violations appear it is because the given nu- 
merical method allows for them. 

Therefore, I argue that it is desirable to control the dis- 
cretized constraints, taking into account possible growth 
introduced both by the formulation of the equations, and 
by the numerical method. The purpose of this paper is 
to present a method for doing so. 

As one is interested in solutions that satisfy the con- 
straints, the Einstein evolution equations can be written 
in infinitely different ways by changing the system's be- 
havior off of the constraint surface. For definiteness I will 
now concentrate on first order (both in space and time) 
formulations, though the ideas of this paper clearly do 
not depend on this and can be applied to second order 
formulations as well. Beginning with any formulation, 
e.g. a symmetric hyperbolic one, one can add to the right 
hand side (RHS) terms proportional to the constraints. 
That is 



E 



{u, t, x)djU -f B(u, t, x) + /iC 



(1) 



where B and u are a vector valued functions (containing 
the metric and related variables), 



C= (Ci...C„) 



(2) 



with C'i ~ C'i{u,dju), is a vector valued constraint (con- 
taining the Hamiltonian, momentum, and perhaps other 
non-physical constraints that appear as integrability con- 
ditions when the system is written in first order form), 
and A, arc matrix-valued functions. The constraints 
(0) in Einstein's equations are quasilinear. That is, non- 
linear in the field variables u, but linear in the spatial 
derivatives, diU. Adding the constraints in a linear way. 
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as in equation yields a quasilinear system of hyper- 
bolic partial differential equations, for which many re- 
sults from the mathematical literature apply [llj. 

The freedom in choosing fi was originally used to cast 
Einstein's equations in symmetric hyperbolic form, later 
to obtain "physical" characteristic speeds (see 3 for re- 
views) , and more recently to improve the stability proper- 
ties of formulations around fixed backgrounds (the num- 
ber of papers in the area is too large to be reviewed here; 
see 0j and references therein). Typically ^ is chosen as 
a constant matrix; however, it does not appear to be the 
most effective way of controlling deviations off the con- 
straint surface in a generic case for two reasons: i) some 
information about the solution (typically, the assumption 
that it will be close -or actually identical - to some known 
background) is usually needed in order to choose an opti- 
mal fj,. However, this kind of information will generically 
not be available, ii) In principle, constants do not seem 
the best way to control varying fields. 

Choosing constant values for ^ seems to be favored for 
historical reasons, and it is actually not required. Sym- 
metric hyperbolicity and physical characteristic speeds, 
for example, can be achieved even if using functions of 
spacetime, provided they are a priori given. See 6] for 
one example. 



II. DYNAMICAL MINIMIZATION OF THE 
CONSTRAINTS GROWTH: THE MAIN IDEA 

Consider the constraints C = C{u, Diu), where the op- 
erator D may be either a continuum or discrete deriva- 
tive. Later D will represent discrete derivatives, since the 
motivation is to control not only continuum constraint 
violations, but also violations that arise in the discrete 
equations. For any slicing of spacetime St x TZ, define a 
norm Nc — Nc{t) for the constraints. Note I will always 
refer to this norm, not the norm associated with the main 
field variables (will comment on this in the last section). 
For simplicity I choose the L2 norm. 



C 



where integration is on the spatial hypersurface St • 
One wants this norm either not to grow during an evo- 
lution, or least to grow as slowly as possible. As C is a 
function of the field variables, any solution to the evolu- 
tion equations automatically determines the rate growth 
for Nc, independent of whether one is considering the 
fully discrete, semidiscrete or continuum systems. I will 
concentrate on the last two cases (i.e., with continuous 
time) . In principle the dependence of Nc on the evolution 
equations is complicated. However, given that /i is added 
linearly to the right hand side of the evolution equations, 
as in EqC] the time derivative of Nc can be written as 



Nc = 



St 



St 



da . 



da 



DkUj 



Using the evolution equations Eq.QJ and allowing fi to 
depend only on time (not on space), the time derivative 
of the norm can then be written as 



iVc = -I- trace(^ x J^^)) 



(3) 



where 



7(1) 



r(2) 



E(^'A 



da 



E 



da. 



dD, 



-Dk 



da 

duj 



a 



(4) 
(5) 



The quantities {iVc, I^^-* , /'■^-', //} depend only on time, in 
the case of /i by assumption and in the others because 
spatial integration is carried out. The splitting is a 
crucial step in the method presented here, since one gains 
analytical control on the norm growth as a function of /i. 
More explicitly: Nc is a linear function of fi (because 
the constraints are added linearly to the evolution equa- 
tions). 

The idea now is to choose /j, such that Nc has some de- 
sired behavior within the freedom allowed by a given for- 
mulation of the evolution equations. As a simple exam- 
ple, consider a single scalar equation. One could choose 
TVc = by defining /i = — J^^^/I^^). Or one could choose 
the norm to decay at a constant or exponential rate by 
defining = -(/(i)+d2)//(2) or = -(/(i)+d2iV,)//(2), 
respectively, with d some constant. On the other hand, 
if the different components of the /i matrix are restricted 
to some sets (this is a restriction that appears quite often 
when requiring symmetric hyperbolicity, physical charac- 
teristic speeds, or both) one could define /i(t) as the one 
that gives the "smallest" growth of Nc{t) in those sets, 
etc. 

This could be done at the continuum. However, the 
problem is that /i would then depend on the constraints, 
which in turn depend on the derivatives of the field vari- 
ables. As a consequence, the resulting evolution equa- 
tions would not be quasilinear, and it is far from clear 
that they would define a well posed initial, or initial- 
boundary, value problem. Indeed, in Section IV I present 
numerical evidence that strongly suggests that the result- 
ing system would be ill posed. 

Therefore, I propose to define /i through a single reso- 
lution run and interpolate the obtained discrete function 
in order to have it defined at all times. In an actual 
computation one would choose a fixed resolution and dy- 
namically compute l'-^-',/^^-' and thus a that gives the 
desired norm behavior, making sure that the resulting 
/i(t) is within the range allowed by symmetric hyper- 
bolicity (symmetric hyperbolicity is not a requirement 
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of the method, but ensures well posedness of the initial- 
boundary value problem, provided appropriate boundary 
conditions are given). In order to ensure numerical sta- 
bility this /i matrix valued function has to be kept fixed 
at other resolutions, /i will depend on the original reso- 
lution used to define it, but this is not a problem. The 
important thing is to keep the constraints under control 
for a given resolution. They will also remain under con- 
trol for better resolutions if one is in the convergence 
regime. 

/i will also depend on the given problem of interest, 
and on the numerical method used to solve the equa- 
tions. That is, for any chosen desired behaviour for Nc, 
there will be one fj., and one well posed formulation of the 
problem for each physical situation one wants to solve for 
with a given numerical method. This is not a practical 
problem and, in fact, it seems difficult to find a formula- 
tion that has optimum stability properties for all possible 
solutions and all possible numerical methods. 



III. CONTROLLING THE CONSTRAINTS 
NORM GROWTH: MORE DETAILS 

If possible, one might want the constraints norm 
growth to be identically zero at all times Nc = 0- A po- 
tential problem with this approach is that the resulting fi 
might have large values, rapid variations in time, or both. 
The resulting set of PDE's would then be stiff, having 
equivalent short timescales, in which case a very small 
Courant factor would be required. In the next section 
I show an example where this happens and argue that 
better results may be obtained by not enforcing Nc = 0. 
Even if one has enough freedom to enforce Nc = 0, it 
appears that allowing Nc for some fluctuations around 
the initial value gives evolutions with better stability. I 
now elaborate on one way of doing so, assuming one has 
enough freedom in fi. This approach is highly non unique, 
and many variations of the fundamental idea presented in 
this paper seem possible, and worth further exploration. 
If one chooses 

Nc = -aNc (6) 

with a > 0, any violation in the constraints will decay 
exponentially 



Nc{t + At) = Ncit) 



-aAt 



(7) 



Choosing a constant a makes the constraints decay for all 
times, I will discuss this in the last section. Here I take 
a different approach. I choose a a tolerance value for the 
norm, Nc = T, and solve for /i such that the constraints 
decay to this tolerance value after a given relaxation time. 
More precisely, I choose a such that after time UaAt the 
constraints will have the value T. Replacing Nc{t + At) 
by T in equation ((TJ and solving for a gives 



ait) 



1 



UaAt 



■In 



T 



I now solve Nc = —aNc = + trace(/x/(^^) for fj,. One 
solution (non unique, since the equation is scalar and fi 
is a matrix) is, assuming J*-^) is invertible, and ommiting 
all the indexes. 



[/(^)]"'(ajV, + /W) 
dim [/(2)] 



(9) 



with a{t) given by eq.®. Thus, for any value of Nc at 
time t, at time t + UaAt the value will be T. Later in this 
paper At will be the discrete time step and, in particular, 
if rta=l, then Nc = T at each time step. 

It must be emphasized that these results hold in the 
semidiscrete case. I.e., for any - not necessarily high 
- spatial resolution. But since time has been assumed 
continuous, the Courant factor must be such that the 
fully discrete simulation faithfully represents the above 
semidiscrete calculations (another option would be to 
perform a fully discrete analysis). In the next section 
I present some numerical experiments to study, among 
other issues, the extent to which the fully discrete sys- 
tem represents the semidiscrete analysis. As we will see, 
at least in the examples here considered, standard values 
of the Courant factor already yield good results. 



IV. SOME NUMERICAL EXAMPLES 

In these numerical experiments I use the method of 
lines and finite differencing. The difference operator cho- 
sen is the simplest one that satisfies summation by parts 
(second order in the interior, and first order in the bound- 
aries), numerical dissipation is introduced taken into ac- 
count modifications at boundaries, and third order Runge 
Kutta is used as time integrator (see 0). 



A. The equations 

I use the standard ADM equations in spherical sym- 
metry with the exact lapse-area locking choices for the 
lapse and shift ^ to illustrate and study the method 
presented here. This choice of lapse and shift amounts 
to specifying in an arbitrary, but a priori, way the lapse 
as a function of spacetime, and the shift as /3 = arKb. 
This choice of shift implies that the radial area does not 

1/2 

change in time, ggg = 0, and one can choose r = ggg as 
a coordinate. The 3-metric and extrinsic curvature in co- 
ordinates r, 9, (p then take the form g.^ = diag{a^, , r^), 
K'^^diag{Ka,Kh,Kb). 

The evolution equations, with the product of the 
Hamiltonian constraint and fi — fi(t) added to the time 
derivative of the Ka equation, are 



Nc{t) 



(8) 



a = ariaKb)' + a [a{Kb - Ka) + rKba] (10) 
= (a'r + 2a)a^^r^'^a - a^^a" + 

a[rKbK + Ka{Ka + 2Kb)]+ fiH (11) 
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a (a 



[Kb{Ka + 2Kk)+r-^{l-a-^)\ (12) 



Constraints norm versus time 
No constraints control 



These equations constitute a strongly hyperbolic sys- 
tem, and the characteristic speeds are "physical" (along 
the light cone or normal to the spatial hypersurfaces): 
/?, /? ± a/a, for any function /i(t). That is, the equations 
are strongly hyperbolic even for /i = 0, i.e., the standard 
ADM equations in spherical symmetry with this choice 
of gauge are strongly hyperbolic. 

The Hamiltonian and momentum constraints are, re- 
spectively, 

H = -X,a' + ^{l-a-'') + 2Kb{Ki, + 2Ka) (13) 

M = K+ ^''' (14) 
r 

I perform these experiments with the Painleve- 
GuUstrand (PG) slicing of the Schwarzschild black hole 
spacetime. 



a = 1 



2r 



a — 1 , (3 — arKb 



(15) 



(16) 



This is a stationary, exact solution, of the evolution 
and constraint equations, (|10I11I12|I . and H13I14|) . respec- 
tively. 

As boundary conditions I set the time derivative of the 
incoming characteristic modes to zero. This will define 
a well posed initial-boundary value problem but in gen- 
eral will violate the constraints. One possibility would 
be to derive constraint-preserving boundary conditions 
for this problem. However, this is not the purpose of this 
paper but, rather, to devise a mechanism to control the 
constraint growth having fixed, among other things, the 
boundary conditions. This is not meant to be a replace- 
ment for constraint-preserving boundary conditions but, 
instead, complementary to them. 



B. Numerically stable but with fast growing errors 
simulations 




FIG. 1: Evolution of the Painlevee-Gullstrand initial data, 
without controlling the constraints (jj. — 0). Nc converges to 
zero when resolution is increased, but at fixed resolution it 
grows fast in time, making the code crash around t — lOAI. 



One could attempt to modify the numerics in order 
to minimize the growth of the constraints, but I do not 
pursue this approach here. 



An interesting feature of the instabilities of the sim- 
ulations shown in Figure 1 is that the norm growth is 
initially negative, but it becomes positive after a single 
time step, triggering the instability, see (|17|l (there, as 
in the rest of the paper, Nc is computed through the 
semidiscrete expression (|2J). 



Ar 


At 


NM 


iVc(Ai) 


M/5 


5.0 X 10^2 


-5.01 X 10-3 


1.47 X 10-2 


M/10 


2.5 X 10^2 


-1.81 X 10-3 


3.85 X 10-3 


M/20 


1.25 X 10-2 


-5.39 X lO-'^ 


9.54 X 10-4 


M/40 


6.25 X 10-3 


-1.46 X IQ-'^ 


2.34 X 10-4 


A//80 


3.12 X 10-3 


-3.82 X 10-5 


5.76 X 10-5 



(17) 



Evolutions of the PG spacetime with equations 
(|10I11I12|1 can give rise to errors that grow fast in time, 
even with a stable numerical method. 

Figure 1 shows Nc vs. time for evolutions of PG initial 
data, at different resolutions. Typical numerical parame- 
ters are chosen: Courant factor A = 0.25, inner and outer 
boundaries at r.i — M, Tq = 40M, respectively, and dis- 
sipation parameter cr = 0.1. The constraints converge to 
zero with resolution, but at fixed resolution they grow 
very fast in time. This is not a peculiarity of this formu- 
lation and this numerical method but, indeed, is typical 
for free evolutions of Einstein's equations in the strong 
field regime. 



A static choice of /i cannot account for these fluctua- 
tions. In particular, choosing /i using only the growth 
rate around the background solution, which in this case 
would be the initial data, and not the errors introduced 
during numerical integration, in general will not be able 
to correct these kind of fluctuations in Nc- 



C. Strictly enforcing semidiscrete norm 
preservation 

I now discuss some numerical experiments using a 
dynamically calculated /i(i) such that TVc = in the 
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Constraints norm versus time 

Different ii used, obtained by conservation of ttie semidiscrete constraints norm 




Convergence test 
H defined by enforcing semidiscrete constraints norm preservation 
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FIG. 2: Strict conservation of the semidiscrete norm {Nc = 0), 
using different resolutions to define i.e. /i is dynamically 
computed at each resolution. The fact that the errors increase 
when resolution is increased strongly suggest that defining fi 
this way leads to an ill-posed problem. The appropriate way 
of doing a convergence test is shown in Figure El 

semidiscrete case. For a chosen resolution, /i is defined 
at each time step by 

f -/(i)(7iAi)//(2)(nAt) , if /(2)(nAt) 7^ 
/i(nAi) ^ I fi{{n- l)At) , if /(2)(nAt) = and n > 1 
[ . if (nAt) = and n = 

(18) 

and iJ,{nAt) is in turn used to compute the RHS needed 
to advance u (nAt). 



Figurc|21shows the resulting norm as a function of time, 
for evolutions enforcing A^c = at different resolutions. 
The code blows up rather fast, and this happens at ear- 
lier times when resolution is increased. This last feature 
should not be seen as an indication of ill-posedness of the 
method, since Figure |21 should not be seen as a conver- 
gence test. For such a test one should dynamically define 
fi(t) with a given resolution and then keep that fi fixed 
for all other runs. Figure El shows the result of doing so, 
defining fj,{t) through a run with Ar = M/5 and then 
keeping it fixed for a convergence test . As expected, 
and in contrast to Figure[21 the norm does decrease as the 
resolution is increased. Figure (21 strongly suggests that 
if one insisted in dynamically defining fi at each resolu- 
tion, instead of fixing it at a single resolution, one would 
end up with an ill-posed problem. On the other hand, 
fixing ^ leads to a strongly hyperbolic problem and no 
convergence problems appear. 



Considering these runs with Nc = 0, it appears that 
enforcing semidiscrete norm preservation is "too rigid" 



FIG. 3: Strict conservation of the semidiscrete norm, defining 
by a Ar = M/5 run and then fixing /i for all other runs. 



|i functions 
Conservation of the semidiscrete constraints norm 




time [M] 



FIG. 4; This figure shows fj,{t) for the runs of Figure 2. 



(recall the discussion at the beginning of Section III). 
Figure 01 shows the functions obtained through the 
runs of Figure [3 Notice that when increasing resolu- 
tion the fi that is needed in order to preserve the norm 
increases quite fast in absolute value. For example, with 
Ar = M/80 the initial values of fi are of order 10^. Since 
fi appears in the principal part of the equations, having 
a large value implies the need of a small Courant factor 
to follow the simulation; while Figures (21 and 01 on the 
other hand, where obtained using the same Courant fac- 
tor when changing resolution. A Fourier decomposition 
of the numerical solution, or an explicit von-Neumann 
analysis of a linear, constant coefficient problem could 
give further insight into this question. 

A second possibility would be to perform a fully dis- 
crete analysis (as opposed to the semidiscrete one used in 
this paper) , and to enforce strict conservation of the fully 
discrete norm, but these issues are not pursued here. 
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10 



H used for convergence test 

Obtained by setting a tolerance value 



Tol = .00051 




10 10 
(tirae+1) [M] 



10 



Convergence test 

defined by setting a tolerance value 




10 



time [M] 



FIG. 5: This figure shows ^ calculated by setting a tolerance 
value, T = 5 X 10~^, and requiring the semidiscrete norm to 
return to this value after every timestep. is defined using 
Ar = M/5, for which the initial, truncation value of the norm 
is Nc = 4.8 X 10~*. 



FIG. 6: A convergence test done with /i shown in Figure |3 
Compare the lifetime of these simulations with those of Figure 
□ 

Norms versus time 

For the constraints and for the main field variables 



D. Dynamically achieving a tolerance value 



Norm for the constraints 

Norm for the main field variables 



I now set a tolerance value T, and choose jjL such that 
Nc = T after a given number of timesteps. In doing so 
I use the semidiscrete expressions m . If IS zero, 
fi will be copied from its previous value, or set to zero if 
this happens in the first iteration. 



Figure [51 shows results with a dynamically defined /i, us- 
ing the same resolution, domain, Courant and dissipation 
factors as before. Ha is chosen to be one, corresponding 
to the norm returning to the specified tolerance value 
after every time step. 

The discretized constraints are, of course, not zero ini- 
tially, even when the initial data is analytic, because 
computing them involves finite differencing the initial 
data. The initial value for the norm for this resolution 
is Nc — 4.8 X 10"''. Therefore I choose a tolerance value 
T = 5 X 10"**, to keep the constraints close to its initial 
value. 

Figure shows the resulting norms for a convergence 
test performed with the same resolutions used in Figure 
n and the fi shown in Figure |31 The code runs for lO^M 
without any sign of instabilities, even with the coarse 
resolution that was used to define /i. 



10 10 
(time+l) [M] 



FIG. 7; Norms associated with the constraints and with the 
main field variables, for the run used to obtain /i shown in 
Figure |^ 



In these experiments, not only does the norm associated 
with the constraints remains close to its initial, trunca- 
tion value, but the same happens with the norm associ- 
ated with the main field variables, N = J (a'^ + K^ + K^). 
Since the exact solution is stationary, N being close to its 
initial value means that the method does not introduce 
new errors. Figure [7| shows Nc{t) and N(t) for the run 
used to define /i. 



1. Dependence of ^ on the tolerance value 



It could happen that one controls the constraints at the 
price of introducing other errors. This does not happen, 
at least in the numerical experiments here considered. 



In the previous examples I chose the tolerance value, 
T, to be roughly the initial, truncation error for Nc- Al- 
though this seems the natural thing to do, I will now 
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Constraints norm versus time 

Same resolution, different tolerance values 



Constraints norm versus time 

Same resolution, different tolerance values 



Tol=0 . 1 
Tol=0 .05 
Tol=0 .01 
Tol=0 .005 
Tol=0 .001 
Tol=0 .0005 
Tol=0 .0001 




time [M] 

FIG. 8: This figure siiows tlie sensitivity of the method to the 
chosen tolerance value T. Values too large (a couple of order 
of magnitudes bigger than the initial truncation error A^c(O)), 
or smaller than Nc{0), lead to instabilities. Otherwise, no fine 
tuning of T is required. 

|i versus time 
Too large and too small tolerance values 



Tol=0 . 0001 
Tol=0 . 1 




(time+1) [M] 



FIG. 9: fi for two unstable cases, with too large {T — 10 ^) 
and too small tolerance values (T = 10"''). 



discuss sensitivity of the stability method with respect 
to the chosen value of T. 



Figure |S1 shows results when choosing different tolerance 
values to define fi, but keeping all other numerical pa- 
rameters (in particular, the resolution, i.e. Figure |S1 is 
not a convergence test) fixed. When T is much larger 
than the initial truncation value of Nc, large errors allow 
triggering of nonlinear instabilities. On the other hand, 
values of T considerably smaller than the initial trunca- 
tion one forces /i to change very fast, and to large values. 
In either case the code crashes. The /i functions for two 
such cases are shown in Figure |3 



le+Ol 



Tol=0.1 
Tol=0.0 
Tol^O.O 
Tol=0.0 
Tol=0.0 
Tol-0.0 
Tol=0.0 




(time+1) [M] 

FIG. 10: The figure shows details of plot|Hl near t = 0. The 
norms in this experiment are the same in the initial data and 
first time step, as explained in the body of the paper. 



Figure ^| shows the details of Figure |S1 near t — 0. 
Notice that not only all the runs begin with the same 
(truncation) value for Nc (this is expected, since the same 
resolution is used), but Nc is also the same in the very 
first timestep. The reason for this is that fj, = initially 
for all the runs, because /^^^ = in the initial data. This 
is a coincidence of the initial data being evolved in this 
example: the only discretization error in the computa- 
tion of the Hamiltonian constraint in this model comes 
from the finite differencing of a, which in this example is 
initially 1 for all grid points. 



E. Intentional constraint violations 

As a final numerical example, now I consider initial 
data that violates the constraints at the continuum. The 
equations still use the exact lapse and area-locking shift 
given by Ea. (|15|l . but now the initial data are given by 



Ka 

Kb 



(l + Sie-('-'-«)') 
f3 



2r 



1 + 526" 



(r-ro)' 



r 



(l + =S3e-('-'^°)') 



(19) 
(20) 
(21) 



When Si = 0, this data reduces to that of the PG black 
hole. The constants Si are used to introduce an "asym- 
metry" in the perturbations of the PG initial data. The 
perturbations are rather large, with si or order unity. 
In this example, si = 1,S2 = 0.5, S3 = 0.7. When 
/I = the code crashes before t — AM at a resolution 
Ar — Mjh. At this resolution Nc for the initial data is 
iVc(O) = 2.66 X 10~^ (one order of magnitude larger than 
the associated truncation error for the case Si = 0, which 
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FIG. 11: The norm Nc, with and without control of the con- 
straints, for highly perturbed initial data, at a resolution of 
Ar = M/5. The initial, truncation error for Nc for this reso- 
lution is Nc{0) = 2.66 X 10~^, and the tolerance value is set 
to T = 5 X 10"^ 
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FIG. 12: n for the plot of Figure lTTl with dynamical controlled 
constraints. After some time fi settles down to the value ^ — 
-5.8 X 10"^ 



for this resolution is Nc = 4.8 x 10^* ). Thus, I choose 
a tolerance value of T = 5 x 10"'^, ni — 10 in order to 
avoid fast variations, and rerun, dynamically obtaining 
Results for the unmodified equations (/i = 0), and 
with a dynamical /i are shown in Figures ^2 ^1 The 
difference in stability is quite striking, especially consid- 
ering that the perturbation of the PG initial data is so 
large. 



V. REMARKS 

There is a lot of interest in finding formulations of Ein- 
stein's equations with good stability properties off of the 
constraint surface. So far, most, if not all, of the work has 
been concentrated in stability properties around known 
solutions. The purpose of this paper has been to intro- 
duce new ideas that hopefully will be helpful in generic 
evolutions, where one has no a priori knowledge of the so- 
lution. This is one fundamental difference with previous 
approaches. 

Another difference is that previous work has concen- 
trated on stability properties at the continuum. Al- 
though this is an important issue, controlling constraint 
violating modes introduced by the numerical method is 
also crucial. Thus, this proposal aims to control discrete 
constraint violations introduced both by the underlying 
formulation of equations at the continuum, and whatever 
numerical method one has chosen. 

Sometimes it seems that one would ideally want the 
constraints to decay exponentially to zero Al- 
though this might be attractive at the continuum, in 
what concerns the discrete constraints it might be bet- 
ter to maintain them around the initial, truncation error. 
For example, beginning with initial data corresponding to 
a stationary spacetime, it seems difficult to imagine that 
by making the constraints decrease in time through nu- 
merical integration (keeping resolution fixed), one would 
end up with a numerical solution that has smaller errors 
than the initial data, which was computed by just numer- 
ically evaluating the exact solution at given gridpoints. 
Thus, I have concentrated here on keeping the discrete 
constraints at some tolerance value, instead of forcing 
them to decay to zero. However, there may be other 
scenarios where it is preferable to make them decay to 
values smaller than the initial one, and the method here 
presented allows for this as well. Along these lines, an ex- 
tension of the work presented here is to use the evolution 
equations applied to initial data off the constraint surface 
to produce solutions of the constraints as an alternative 
to relaxing the constraints themselves )14] . 

In some other approaches , a norm for the main vari- 
ables, say N{t) = Jg^J^i'^h is minimized with respect 
to constraint violating perturbations (around a given, 
known background) . I could have here chosen to dynami- 
cally minimize this quantity. One potential problem with 
this is that it is not clear how much minimization makes 
the solution closer to the "exact" one. For example, in 
the spherically symmetric example here considered the 
time derivative of the norm for the main variables has 
the form 

TV = ja) + ^/r2) 

Since /i is completely free in this model, one is able to 
make N, for some chosen but arbitrary resolution, as 
negative as one wants. What happens then is that all the 
field variables decay to zero (numerical experiments that 
I have done do confirm this). Clearly this is not what 
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one wants, since the "exact" solution is not identically 
zero. On the other hand, one could also attempt to keep 
N close to its initial, truncation value. However, this 
would be a good idea only if the solution is stationary, 
since otherwise there is no reason for the exact N{t) to be 
close to A^(0) for all times. Thus I have here considered 
a norm associated with the constraints, since they are 
analytically zero, regardless of the solution. 
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